Validation with mortality rates in NFI plots

Introduction

In this report, we aim to estimate the association between the genomic offset (GO) predictions and mortality rates in natural populations of maritime pine, used as a proxy of the absolute fitness of the populations. A positive association between GO predictions and the mortality rates would suggest that GO approaches capture maladaptation in the face of demographic complexity (i.e., the effects of other processes on spatial variance in allele frequencies, such as expansion history and gene flow) as well as the genetic architecture of climate adaptation (polygenic trait architectures, G×E interactions, nonadditive genetic variance).

For that, we use mortality data from the French and Spanish National Forest Inventories (NFI) harmonized in Changenet et al. (2021). The French data relies on temporary plots sampled between 2005 and 2014 while the Spanish data relies on permanent plots sampled during the second (from 1986 to 1996) and third NFIs (from 1997 to 2008). A tree was recorded as dead if its death was dated at less than 5 years ago in the French NFI or if it was alive in the second inventory but dead in the third one in the Spanish NFI.

The different census intervals among the plots as accounting for by modelling the proportion \(p_{i}\) of maritime pines that died in the plot \(i\) during the census interval \(\Delta_i\) with the complementary log-log link and an offset on the logarithm of the census interval \(\Delta_{i}\) for the plot \(i\). The modeling approach was evaluated with simulations in the report 14_ValidationNFI_ModelComparison.qmd.

In the present report, we evaluate the relationship between mortality rates in the NFI plots and the G0 predictions with the two following models.

Model with competition as covariate

This model corresponds to the model M3 in the report 14_ValidationNFI_ModelComparison.html.

\[\begin{align*} m_i &\sim \text{Binomial}(N_i,p_i)\\ \text{log}(-\text{log}(1-p_i)) &= \beta_{0,c} + \beta_{C} C_i + \beta_{GO} GO_i + \text{log}(\Delta_i) \\ \end{align*}\]

• \(N_{i}\) is the total number of maritime pines in the plot \(i\),

• \(m_{i}\) is the number of maritime pines that died during the census interval \(\Delta_{i}\) in the plot \(i\),

• \(\beta_{0,c}\) are country-specific intercepts that account for the methodological differences between the French and Spanish inventories that may bias the estimations,

• \(C_{i}\) is the basal area of all tree species confounded in the plot \(i\) (to account for the competition among trees), with its associated coefficient \(\beta_{C}\),

• \(GO_{i}\) is the genomic offset predicted in the plot \(i\), with its associated coefficient \(\beta_{GO}\).

Model with competition and tree age as covariates

This model corresponds to the model M6_2 in the report 14_ValidationNFI_ModelComparison.html.

\[\begin{align*} m_i &\sim \text{Binomial}(N_i,p_i)\\ \text{log}(-\text{log}(1-p_i)) &= \beta_{0,c} + \beta_{C} C_i + \beta_{DBH} DBH_i + \beta_{C*DBH} C_i DBH_i + \beta_{GO} GO_i + \text{log}(\Delta_i) \\ \end{align*}\]

• \(DBH_{i}\) the mean diameter at breast height (DBH) of the maritime pines in the plot \(i\) (including adults, dead trees and recruitment trees) and its associated coefficient \(\beta_{DBH}\). This variable is expected to account for the different mean tree ages across plots that may impact mortality rates.

We present the results from this model in the main manuscript.

Priors

We used the following weakly informative priors:

\[\begin{bmatrix} \beta_{0,c} \\ \beta_C \\ \beta_{GO} \\ \beta_{DBH} \\ \beta_{C*DBH} \end{bmatrix} \sim \text{Normal}(0,1)\]

NFI data

NFI data was harmonized in Changenet et al. (2021).

Code

nfi_data <- readRDS(here::here("data/NFIdata/NFIrawdata_Changenet2021.rds")) %>% 
  droplevels() %>% 
  as_tibble() %>% 
  dplyr::select(plotcode,
                longitude,
                latitude,
                country,
                yearsbetweensurveys, # number of years between surveys
                surveydate1, # date first survey Spanish NFI
                surveydate, # year of the second survey in the Spanish NFI and the unique survey in the French NFI
                treeNbrJ.M, # number of dead trees in the plot
                treeNbrJ.IMall, # total number of trees in the plot
                BA.ha.plot.1, # basal area of all tree species in the plot
                dbhJ.IMall.plot.mean # mean DBH of the maritime pines in the plot including dead trees but also saplings (recruitment) 
                ) %>% 
  dplyr::rename(nb_years = yearsbetweensurveys,
                nb_dead = treeNbrJ.M,
                nb_tot = treeNbrJ.IMall,
                basal_area = BA.ha.plot.1,
                mean_DBH = dbhJ.IMall.plot.mean,
                first_survey = surveydate1,
                second_survey = surveydate) %>% 
  dplyr::mutate(first_survey = case_when(!is.na(first_survey) ~ str_sub(first_survey,1,4) %>% as.numeric()), # keep only the year of the 1st survey
                prop_dead = nb_dead/nb_tot,  # proportion of dead trees
                annual_prop_dead = (nb_dead/nb_tot)/nb_years) # annual proportion of dead trees

Variables in the dataset

• Plot code: plotcode
• Longitude and latitude of the plot: longitude and latitude.
• Country in which the plot is: country (ES = Spain, FR = France).
• The number of years between surveys in the Spanish inventory (which is equal to 5 in the French inventory as mortality is estimated in the five years before the survey date): nb_years.
• The dates of the first and second survey in the Spanish inventory: surveydate1 and surveydate2.
• The year of the second survey in the Spanish inventory and of the unique survey in the French inventory: surveydate.
• The number of dead trees in the plot: nb_dead.
• The total number of trees in the plot: nb_tot.
• The basal area of all tree species in the plot (proxy of the competition among trees): basal_area.
• The mean DBH of the maritime pines in the plot including dead trees but also saplings (recruitment): mean_DBH.

12610 plots in this dataset: 4097 from the French NFI and 8513 from the Spanish inventory.

Filtering

How many NAs in each column?

Code

nfi_data %>% 
  summarise(across(everything(), ~ sum(is.na(.)))) %>% 
  pivot_longer(everything(), names_to = "Variable", values_to = "Number of NAs") %>% 
  kable_mydf(boldfirstcolumn = T)

Variable	Number of NAs
plotcode	0
longitude	0
latitude	0
country	0
nb_years	0
first_survey	4097
second_survey	0
nb_dead	0
nb_tot	0
basal_area	277
mean_DBH	664
prop_dead	664
annual_prop_dead	664

We extract the plots with missing data for the proportion of dead trees and we look at the first eight plots.

Code

df_na_prop_dead <- nfi_data %>% 
  dplyr::filter(is.na(prop_dead))
df_na_prop_dead[1:8,]  %>% kable_mydf(boldfirstcolumn = F)

plotcode	longitude	latitude	country	nb_years	first_survey	second_survey	nb_dead	nb_tot	basal_area	mean_DBH	prop_dead	annual_prop_dead
ES100047A1	-6.34	40.42	ES	11	1990	2001	0	0	2.40	NA	NaN	NaN
ES100094A1	-6.22	40.39	ES	11	1990	2001	0	0	3.52	NA	NaN	NaN
ES100129A1	-6.28	40.37	ES	11	1990	2001	0	0	NA	NA	NaN	NaN
ES100151A1	-6.30	40.37	ES	11	1990	2001	0	0	NA	NA	NaN	NaN
ES100203A1	-6.27	40.35	ES	11	1990	2001	0	0	6.11	NA	NaN	NaN
ES100206A1	-6.52	40.34	ES	11	1990	2001	0	0	NA	NA	NaN	NaN
ES100232A1	-6.17	40.44	ES	11	1990	2001	0	0	NA	NA	NaN	NaN
ES100307A1	-6.93	40.25	ES	11	1990	2001	0	0	NA	NA	NaN	NaN

No maritime pine was recorded in these eight plots (nb_dead = 0 and nb_tot = 0). We check that this is the case for the 664 plots by looking at the unique values in the columns nb_dead and nb_tot in the subset of plots with missing data for the proportion of dead trees (the df_na_prop_dead dataset).

Code

df_na_prop_dead %>% 
  summarise(across(nb_dead:nb_tot, ~ unique(.))) %>% 
  pivot_longer(everything(), names_to = "Variable", values_to = "Unique values") %>% 
  kable_mydf(boldfirstcolumn = T)

Variable	Unique values
nb_dead	0
nb_tot	0

We remove these 664 plots.

Code

nfi_data <- nfi_data %>% 
  dplyr::filter(!is.na(prop_dead))

We check that there are no more plots without maritime pine.

Code

nfi_data %>% dplyr::filter(nb_tot==0) %>% nrow()

[1] 0

Missing data for basal area:

Code

nfi_data %>% 
  summarise('Number of NAs for basal area' = sum(is.na(basal_area))) %>% 
  kable_mydf(boldfirstcolumn = F)

Number of NAs for basal area
29

Since there are not many plots with missing basal area data, we delete them.

Code

nfi_data <- nfi_data %>% 
  dplyr::filter(!is.na(basal_area))

After these filtering steps, there are 11917 plots left.

We save the file for model comparison and we save a copy for the DRYAD repository.

Code

nfi_data %>% saveRDS(file=here("data/NFIdata/NFIdata_cleaned.rds"))

nfi_data %>% write_csv(here("data/DryadRepo/NFIdata_cleaned.csv"))

Inventory dates

Code

country_names <- c(`ES` = "SPAIN",`FR` = "FRANCE")

ggplot(nfi_data, aes(x=second_survey, fill=as.factor(country))) + 
  geom_histogram(alpha=0.5, position="identity") + 
  theme_bw()  +
  facet_wrap(~as.factor(country), labeller = as_labeller(country_names)) + 
  labs(y="Count of plots",x="Inventory date") +
  theme(legend.position = "none")

Code

ggplot(nfi_data, aes(x=first_survey)) + 
  geom_histogram(alpha=0.5, position="identity") + 
  theme_bw()

Code

nfi_data %>% 
  ggplot(aes(x=nb_years, fill=country)) + 
  geom_histogram(alpha=0.5, position="identity") + 
  labs(y="Count of plots",x="Number of years between surveys") +
  theme_bw()

Climatic data

Generating file for ClimateDT

Code

# NFI plots coordinates
# ---------------------
nfi_climateDT <- nfi_data %>% 
  dplyr::select(plotcode, longitude, latitude) %>% 
  write_csv(here::here("data/nfi_coordinates.csv"))

We are going to extract the climatic data between 10 years before the first survey until the date of the second survey. For the French NFI, as there is only one survey, we will extract climatic data for the 15 years before the survey.

Code

nfi_data %>%
  group_by(country) %>%
  mutate(min_year = case_when(country=="FR" ~ second_survey-15,
                              country=="ES" ~ first_survey-10)) %>% 
  summarise(min_year=min(min_year),
            max_year=max(second_survey)) %>% 
  kable_mydf()

country	min_year	max_year
ES	1972	2008
FR	1990	2014

So, we need annual climatic data at the location of the NFI plots for the period 1972-2014.

Calculate climate for each plot

Once we have the annual climatic data at the location of the NFI plots, we generate two datasets:

• one with the past climate data of the NFI plots, across the period 1901-1950.

• one with the climatic data across a period specific to each plot, corresponding to:

  • The period from 5 years before the first inventory date to the second inventory date for the Spanish NFI.

  • The period from 10 years before the inventory date and the inventory date for the French NFI. A tree is considered dead in the French NFI if its death is estimated as less than 5 years before the inventory date.

We included the 5 years preceding the estimated date of tree death, as the climate of the previous years may have lag effects on tree death.

We load the annual climatic data, which corresponds to point estimate at the location of the NFI plots:

Code

clim <- read_csv(here("data/ClimaticData/NFIplots/ClimateDT_cmip6_GFDL-ESM4_NFIcoordinates_ADJ.csv"), show_col_types = FALSE) %>% 
  dplyr::rename(plotcode=ID,
                longitude=Longitude,
                latitude=Latitude,
                elevation=Elevation,
                year=Year)

We calculate the average climatic conditions across the period 1901-1950:

Code

# Function to calculate the average of the climatic variables at the location of the NFI plots
source(here("scripts/functions/calc_avg_clim_var.R"))

clim_ref <- clim %>% calc_avg_clim_var(ref_period = c(1901,1950), id_spatial_points = "plotcode")

Code

# !!Coding warning!! 
######################

# The coordinates from ClimateDT are not exactly the same as the original ones. 
# Therefore, merging the climatic data with the original data has to be done with the `plotcode` column 
# and not the latitude and longitude of the populations.

# Comparing original coordinates with the ones from ClimateDT
comp_coord <- clim_ref %>% 
  dplyr::select(plotcode,latitude,longitude) %>% 
  dplyr::rename(latitude_climDT=latitude, longitude_climDT=longitude)
comp_coord <- nfi_data %>% 
  dplyr::select(plotcode,latitude,longitude) %>% 
  inner_join(comp_coord, by="plotcode") %>% 
  mutate(diff_latitude=latitude_climDT - latitude,
         diff_longitude=longitude_climDT - longitude) %>% 
  dplyr::filter(diff_latitude != 0 | diff_longitude != 0)

# There are very very small differences
range(comp_coord$diff_latitude)
range(comp_coord$diff_longitude)

Code

# checking time periods for Maurizio
# ==================================

time_periods <- clim_survey %>% 
  dplyr::select(min_year_clim,second_survey) %>% 
  group_by(min_year_clim,second_survey) %>% 
  summarise(count=n()) 

time_periods %>% 
  print(n=nrow(.)) %>% 
  kable_mydf()

time_periods %>% filter(count <10)

time_periods %>% 
  write_csv(here("data/ClimaticData/NFIplots/NFIplots_TimePeriods.csv"))

We then calculate the average climatic conditions across the survey periods specific to each NFI plot:

Code

# We calculate the mean climatic values between:
  # - five years before the first inventory date
  # - the second inventory date
clim_survey <- nfi_data %>% 
  mutate(min_year_clim = case_when(country=="FR" ~ second_survey-10,
                                   country=="ES" ~ first_survey-5))

# we calculate the average climatic conditions for each specific period
clim_survey <- clim_survey %>% 
  group_by(min_year_clim,second_survey) %>% 
  group_split() %>% 
  purrr::map(\(x){
  
min_year <- unique(x$min_year_clim)
max_year <- unique(x$second_survey)

clim %>% 
  right_join(x[,"plotcode"],by=c("plotcode")) %>%
  calc_avg_clim_var(ref_period = c(min_year,max_year), id_spatial_points = "plotcode")
  
}) %>% list_rbind() 


# we save in a list the average climatic conditions across the two periods (reference period and period of the NFI surveys)
list(clim_ref = clim_ref,
     clim_survey = clim_survey) %>% 
  saveRDS(here("data/ClimaticData/NFIplots/NFIclimate.rds"))

We save the averaged climatic data at the locations of the NFI plots in the DRYAD repository:

• dataset ClimateDT_NFIPlots_PastClimates.csv for the averaged climates over the reference period 1901-1950.

• dataset ClimateDT_NFIPlots_SurveyClimates.csv" for the averaged climates over the survey periods specific to each plot.

Code

clim_ref %>% write_csv(here("data/DryadRepo/ClimateDT_NFIPlots_PastClimates.csv"))
clim_survey %>% write_csv(here("data/DryadRepo/ClimateDT_NFIPlots_SurveyClimates.csv"))

Climatic space covered

We look at how the 34 sampled populations cover the climatic space of the NFI plots.

Code

# Selected climatic variables
clim_var <- readRDS(here("data/ClimaticData/NamesSelectedVariables.rds"))

# Reference climate of the NFI plots
clim_nfi <- clim_ref %>% 
  dplyr::select(plotcode,contains("ude"),elevation, all_of(clim_var)) %>% 
  dplyr::mutate(subdf="nfi") %>% 
  dplyr::rename(id=plotcode)

# Gene pools of the populations
gps <- readRDS(here("data/GenomicData/MainGenePoolInformation.rds"))[[1]] %>% arrange(pop)

# Reference climate of the populations
clim_pop <- readRDS(here("data/ClimaticData/MaritimePinePops/ClimatePopulationLocationPointEstimates_ReferencePeriods_noADJ.rds"))[[1]]$ref_means %>% 
  dplyr::select(pop,contains("ude"),elevation, all_of(clim_var)) %>% 
  left_join(gps,by="pop") %>% 
  dplyr::rename(id=pop,
                subdf=main_gp_pop) 

# Combining the dataset of the populations and the NFI plots
df<- bind_rows(clim_pop,clim_nfi) %>% 
  mutate(sudf=factor(subdf,levels=c("nfi",unique(clim_pop$subdf))))


# Mean annual temperature vs annual precipitation
# ===============================================
p <- df %>%
  ggplot(aes(x=bio1,y=bio12,label=id)) +
  geom_point(data=df %>% filter(subdf=="nfi"), color="gray85") +
  geom_point(data=df %>% filter(!subdf=="nfi"), aes(color=subdf),size=4) +
  scale_color_manual(values=c("darkorchid3",
                              "gold2",
                              "navyblue",
                              "turquoise2",
                              "orangered3",
                              "green3")) +
  xlab("Mean annual temperature (bio1, °C)") +
  ylab("Annual precipitation (bio12, mm)") +
  geom_text_repel(data=df %>% filter(!subdf=="nfi")) +
  theme_bw() +
  theme(legend.title = element_blank(),
        legend.position = c(0.82,0.87),
        legend.text = element_text(size=14),
        legend.background = element_blank(),
        axis.text = element_text(size=12),
        axis.title = element_text(size=18))

p %>% ggsave(filename = here("figs/ValidationNFI/ClimaticCoverage_Bio1vsBio12.pdf"),
             device="pdf",
             width = 8,
             height=8)

p

Code

# Precipitation seasonality vs Isothermality    
# ==========================================
p <- df %>%
  ggplot(aes(x=bio15,y=bio3,label=id)) +
  geom_point(data=df %>% filter(subdf=="nfi"), color="gray85") +
  geom_point(data=df %>% filter(!subdf=="nfi"), aes(color=subdf),size=4) +
  scale_color_manual(values=c("darkorchid3",
                              "gold2",
                              "navyblue",
                              "turquoise2",
                              "orangered3",
                              "green3")) +
  xlab("Precipitation seasonality (bio15, coefficient of variation)") +
  ylab("Isothermality (bio3, index)") +
  geom_text_repel(data=df %>% filter(!subdf=="nfi")) +
  theme_bw() +
  theme(legend.title = element_blank(),
        legend.position = c(0.17,0.15),
        legend.text = element_text(size=14),
        legend.background = element_blank(),
        axis.text = element_text(size=12),
        axis.title = element_text(size=18))

p %>% ggsave(filename = here("figs/ValidationNFI/ClimaticCoverage_Bio3vsBio15.pdf"),
             device="pdf",
             width = 8,
             height=8)

p

Code

# Temperature seasonality and Summer heat moisture index    
# ======================================================
p <- df %>%
  ggplot(aes(x=bio4,y=SHM,label=id)) +
  geom_point(data=df %>% filter(subdf=="nfi"), color="gray85") +
  geom_point(data=df %>% filter(!subdf=="nfi"), aes(color=subdf),size=4) +
  scale_color_manual(values=c("darkorchid3",
                              "gold2",
                              "navyblue",
                              "turquoise2",
                              "orangered3",
                              "green3")) +
  xlab("Temperature seasonality (bio4, standard deviation ×100)") +
  ylab("Summer heat moisture index (SHM, °C/mm)") +
  geom_text_repel(data=df %>% filter(!subdf=="nfi")) +
  theme_bw() +
  theme(legend.title = element_blank(),
        legend.position = c(0.16,0.89),
        legend.text = element_text(size=14),
        legend.background = element_blank(),
        axis.text = element_text(size=12),
        axis.title = element_text(size=18))

p %>% ggsave(filename = here("figs/ValidationNFI/ClimaticCoverage_Bio4vsSHM.pdf"),
             device="pdf",
             width = 8,
             height=8)

p

Genomic offset predictions

We load the genomic offset predictions estimated from the different GEAs.

Code

df <- lapply(c("control","cand"), function(x){

list_snps <- list()

list_snps[[x]]$GDM <- readRDS(file=here("outputs/GDM/go_predictions.rds"))[[x]][["go_nfi"]]
list_snps[[x]]$GF <- readRDS(file=here("outputs/GF/go_predictions.rds"))[[x]][["go_nfi"]]
list_snps[[x]]$RDA <- readRDS(file=here("outputs/RDA/go_predictions.rds"))[[x]][["go_nfi"]]
list_snps[[x]]$LFMM <- readRDS(file=here("outputs/LFMM/go_predictions_snpsets.rds"))[[x]][["go_nfi"]]


list_snps <- list_snps[[x]] %>% 
  bind_rows() %>% 
  mutate(method_input = x,
         plotcode = unique(nfi_data$plotcode))

return(list_snps)
}) %>% 
  bind_rows() %>% 
  pivot_longer(cols=c("GDM","GF","RDA","LFMM"), names_to = "method_type",values_to ="GO") %>% 
  bind_rows(tibble(GO = readRDS(file=here("outputs/LFMM/go_predictions_allsnps.rds"))[["go_nfi"]],
                   method_type = "LFMM",
                   method_input = "all",
                   plotcode = unique(nfi_data$plotcode))) %>% 
  mutate(method = paste0(method_type, " (",method_input," SNPs)"))

Relationship btw GO predictions and mortality rates

Code

# Colors for the different sets of GO method (GDM, GF, LFMM and RDA) / SNP set (control and candidate SNPs, and all SNPs for LFMM)
colors_coeff <- c("#004586FF","#FF7F0FFF","#FFD320FF","#579D1CFF","#7E0021FF","#83CAFFFF","#87D180FF","#FEB5A2FF","#F02720FF")

Model with competition as covariate

Stan code

Stan code of the model (model M3 in the report 14_ValidationNFI_ModelComparison.html):

Code

stancode = stan_model(here("scripts/StanModels/ValidationNFI_m3.stan"))
print(stancode)

S4 class stanmodel 'anon_model' coded as follows:
data {
  int N;
  vector[N] log_nb_years;                                                           // Offset to account for different census intervals
  vector[N] C;                                                                      // Proxy of the competition among trees (i.e. basal area)
  vector[N] GO;                                                                     // Genomic offset
  int nb_dead[N];                                                                   // Number of dead trees in the plot
  int nb_tot[N];                                                                    // Total number of trees in the plot
  int<lower=0> nb_country;                                                          // Number of countries
  int<lower=0, upper=nb_country> country[N];                                        // Countries
}
parameters {
  vector[nb_country] alpha_country;
  real beta_GO;
  real beta_C;
}
model {
  nb_dead ~ binomial(nb_tot,inv_cloglog(alpha_country[country] + beta_GO * GO + beta_C * C  + log_nb_years));

  alpha_country ~ normal(0, 1);
  beta_GO ~ normal(0, 1);
  beta_C ~ normal(0, 1);
} 

Running the model

Running the model for each method (GDM, LFMM, GF and RDA) and each SNP set (all SNPs - only for LFMM-, candidate SNPs and control SNPs):

Code

# Parameters to estimate
params_to_estimate <- c("alpha_country","beta_GO","beta_C")

# Running one model for each of the nine sets of method / SNP set
coefftab <- lapply(unique(df$method), function(method_i){
  
# Subset the data - keeping only one set of method / SNP set
    sub_data <- df %>% 
      filter(method == method_i) %>% 
      inner_join(nfi_data, by="plotcode")
      
# Data in a list for Stan 
    stanlist <- list(N = nrow(sub_data),
                     nb_dead = sub_data$nb_dead,
                     nb_tot=sub_data$nb_tot,
                     GO = (sub_data$GO - mean(sub_data$GO))/sd(sub_data$GO),
                     C = (sub_data$basal_area - mean(sub_data$basal_area))/sd(sub_data$basal_area),
                     nb_country=length(unique(sub_data$country)),
                     country=as.numeric(as.factor(sub_data$country)),
                     log_nb_years=log(sub_data$nb_years))

       

# Running the Stan model
mod <- sampling(stancode, data = stanlist, iter = 2000, chains = 4, cores = 4,init=0, save_warmup = FALSE)

# Model coefficients
conf95 <- broom.mixed::tidyMCMC(mod,
                                pars=params_to_estimate,
                                droppars = NULL, 
                                robust = FALSE, # give mean and standard deviation
                                ess = T, # effective sample size estimates
                                rhat = T, # Rhat estimates
                                conf.int = T, conf.level = 0.95 # include 95% credible intervals
                                ) %>% 
    dplyr::rename(conf95_low=conf.low,
                  conf95_high=conf.high,
                  mean=estimate,
                  std_deviation=std.error)
  
  
  broom.mixed::tidyMCMC(mod,
                        pars=params_to_estimate,
                        droppars = NULL, 
                        robust = TRUE, # give median and mean absolute deviation (= avg distance btw each data point and the mean)
                        ess = F, rhat = F, 
                        conf.int = T, conf.level = 0.80 # include 80% credible intervals
                        ) %>% 
    dplyr::rename(conf80_low=conf.low,
                  conf80_high=conf.high,
                  median=estimate,
                  mean_absolute_deviation=std.error) %>%
    inner_join(conf95,by=c("term")) %>% 
    mutate(method=method_i,
           method_type=unique(sub_data$method_type),
           method_input=unique(sub_data$method_input),
           .before=1)
    
  }) %>% bind_rows()

coefftab %>% saveRDS(file=here("outputs/ValidationNFI/ModelM3_coefftab.rds"))

Model coefficients

Below, we plot the 95% credible intervals of:

• the \(\beta_{GO}\) coefficients, which stand for the association between mortality rates in NFI plots and the genomic offset predictions (i.e. capturing the potential maladaptation of the populations at the location of the NFI plots).

• the \(\beta_C\) coefficients, which stands for the association between the tree basal area (all species confounded) in the NFI plots and the mortality rates.

Code

coefftab <- readRDS(file=here("outputs/ValidationNFI/ModelM3_coefftab.rds"))
  
coeff_match <- list(beta_GO=list("Regression coefficients $\\beta_{GO}$",c(0.77,0.85)),#"$\\beta_{GO}$ estimates",
                    beta_C=list("Regression coefficients $\\beta_C$",c(0.77,0.16)))

p <- lapply(c("beta_GO","beta_C"), function(coeff){
  
p <- coefftab %>% 
  filter(term==coeff) %>% 
  mutate(method_input=case_when(method_input=="all" ~ "All SNPs",
                                method_input=="cand" ~ "Candidate SNPs",
                                method_input=="control" ~ "Control SNPs")) %>% 
  mutate(method_input=factor(method_input, levels=c("All SNPs",
                                                    "Candidate SNPs",
                                                    "Control SNPs"))) %>% 
    ggplot(aes(x = method_type,
               y = median,
               ymin = conf95_low, 
               ymax = conf95_high,
               color=method_input,
               shape=method_input)) +
  
  {if(coeff=="beta_GO")  geom_rect(inherit.aes=FALSE, aes(xmin=-Inf, xmax=Inf, ymin=0, ymax=Inf), color="transparent", fill="green", alpha=0.01)} + 
  geom_hline(yintercept = 0,color="gray") +
  geom_pointinterval(position = position_dodge(width = .4),point_size=3.5,size=3) +
  ylab(TeX(coeff_match[[coeff]][[1]])) + xlab("") +
  ylim(c(-0.15,0.23)) +
  scale_color_manual(values=colors_coeff) +
  scale_shape_manual(values = c(rep(16,6),rep(17,3))) +
  theme_bw() +
  labs(color="SNP set",shape="SNP set") +
  theme(axis.text.x =  element_text(size=13),
        axis.text.y = element_text(size=13),
        axis.title.y = element_text(size=16),
        axis.title.x = element_text(size=1),
        legend.title=element_text(size=13), 
        legend.box.background = element_rect(colour = "gray", linewidth=0.6),
        strip.text.x = element_text(size = 16),
        strip.background = element_blank(),
        legend.position = coeff_match[[coeff]][[2]],
        legend.text=element_text(size=12),
        panel.grid.minor.x=element_blank(),
        panel.grid.major.x=element_blank()) +
  guides(color=guide_legend(ncol=1),
         shape = guide_legend(override.aes = list(size =2 )))

# save in pdf and png
ggsave(p,file=here(paste0("figs/ValidationNFI/ModelM3_",coeff,"_IntervalPlots.pdf")),
       device="pdf",
       height=5,
       width=5)

return(p)

})

p

Interpretation

As expected, the basal area is positively associated with mortality rates.

For the four methods (GDM, LFMM, GF and RDA), populations with higher genomic offset based on the candidate SNPs (and all SNPs for LFMM) had higher mortality rates in the NFI plots. Interestingly, this was not the case with genomic offset predictions based on the control SNPs: for GF and LFMM, populations with higher genomic offset based on the control SNPs also showed higher mortality rates but for GDM and RDA, populations with higher genomic offset showed lower mortality rates in the NFI plots. This confirms the importance of selecting a set of alleles potentially under selection before using the genomic offset predictions to predict maladaptation in natural populations (or using all SNPs with the LFMM approach).

We can also note the very high association between mortality rates and the genomic offset predictions based on candidate SNPs and from the GF and GDM methods.

\(\beta_{GO}\) estimates against a null distribution

We evaluated the probability that the estimated relationship between the mortality rates and GO predictions could be obtained by chance. For that, we ran again the model on the NFI data but we replaced GO predictions by a randomly generated variable \(X\) such as \(X \sim \mathcal{N}(0,1)\). We repeated this step 100 times and obtained a distribution of \(\beta_X\) coefficients capturing the relationship between mortality rates and a randomly generated variable. We then compared this distribution of \(\beta_X\) coefficients to the estimated \(\beta_{GO}\) coefficients capturing the relationship between the mortality rates and the GO predictions.

Code

# Load the coefficient values of models with predicted GO
coefftab <- readRDS(file=here("outputs/ValidationNFI/ModelM3_coefftab.rds")) %>% 
  mutate(run_ID=method,
         GO="predicted") %>% 
  dplyr::filter(term == "beta_GO") %>%
  dplyr::select(run_ID, GO, mean, contains("conf"))

# Load the coefficient values of models with random GO and merge with the predicted GO coefficient table
coefftab <- readRDS(file=here(paste0("outputs/ValidationNFI/Simulations/m3_truedata_randomGO_100simulations.rds"))) %>% 
  bind_rows(.id="sim_ID") %>% 
  dplyr::filter(term == "beta_GO") %>%
  mutate(run_ID = paste0("sim ",sim_ID),#run_ID = paste0("Simulation ",sim_ID),
         GO="random") %>% 
  dplyr::select(run_ID, GO, mean, contains("conf")) %>% 
  bind_rows(coefftab)
  
# Plot the mean and 95% credible intervals with interval plots
p <- coefftab %>% 
  ggplot(aes(mean, reorder(run_ID, mean))) +
  geom_vline(xintercept = 0, color="gray30") +
  geom_errorbar(aes(xmin = conf95_low, xmax = conf95_high, color=GO)) +
  scale_color_manual(values=c("magenta","gray"),
                     name = "Genomic offset") +
  geom_point() +
  theme_bw() +
  labs(x="GO effect size", y="") +
  theme(axis.text.y = element_text(size=6))

# save the plot in pdf
ggsave(p,
       file=here(paste0("figs/ValidationNFI/ModelM3_IntervalPlots_RandomGO_vs_PredictedGO.pdf")),
       device="pdf",
       height=10,
       width=8)

p

This graph shows that it is very unlikely that the relationship between the mortality rates and GO predictions of GDM + candidate SNPs, GF + candidate SNPs, RDA + candidate SNPs and GDM + control SNPs are obtained by chance.

Model with competition and tree age as covariates

Stan code

Stan code of the model (model M6_2 in the report 14_ValidationNFI_ModelComparison.html):

Code

stancode = stan_model(here("scripts/StanModels/ValidationNFI_m6_2.stan"))
print(stancode)

S4 class stanmodel 'anon_model' coded as follows:
data {
  int N;
  vector[N] log_nb_years;                                                           // Offset to account for different census intervals
  vector[N] C;                                                                      // Proxy of the competition among trees (i.e. basal area)
  vector[N] DBH;                                                                    // Proxy of the average tree age (i.e. mean DBH)
  vector[N] GO;                                                                     // Genomic offset
  int nb_dead[N];                                                                   // Number of dead trees in the plot
  int nb_tot[N];                                                                    // Total number of trees in the plot
  int<lower=0> nb_country;                                                          // Number of countries
  int<lower=0, upper=nb_country> country[N];                                        // Countries
}
parameters {
  vector[nb_country] alpha_country;
  real beta_GO;
  real beta_C;
  real beta_DBH;
  real beta_C_DBH;
}
model {
  nb_dead ~ binomial(nb_tot,inv_cloglog(alpha_country[country] + beta_GO * GO + beta_C * C  + beta_DBH * DBH + beta_C_DBH * C .* DBH + log_nb_years));

  alpha_country ~ normal(0, 1);
  beta_GO ~ normal(0, 1);
  beta_C ~ normal(0, 1);
  beta_DBH ~ normal(0, 1);
  beta_C_DBH ~ normal(0, 1);
} 

Running the model

Running the model for each method (GDM, LFMM, GF and RDA) and each input (all SNPs - only for LFMM-, candidate SNPs and control SNPs):

Code

# Parameters to estimate
params_to_estimate <- c("alpha_country","beta_GO","beta_C", "beta_DBH","beta_C_DBH")

# Running one model for each of the nine sets of method / SNP set
coefftab <- lapply(unique(df$method), function(method_i){
    
# Subset the data - keeping only one set of method / SNP set
      sub_data <- df %>% 
      filter(method == method_i) %>% 
      inner_join(nfi_data, by="plotcode")
      
# Data in a list for Stan 
    stanlist <- list(N = nrow(sub_data),
                     nb_dead = sub_data$nb_dead,
                     nb_tot=sub_data$nb_tot,
                     GO = (sub_data$GO - mean(sub_data$GO))/sd(sub_data$GO),
                     C = (sub_data$basal_area - mean(sub_data$basal_area))/sd(sub_data$basal_area),
                     DBH = (sub_data$mean_DBH - mean(sub_data$mean_DBH))/sd(sub_data$mean_DBH),
                     nb_country=length(unique(sub_data$country)),
                     country=as.numeric(as.factor(sub_data$country)),
                     log_nb_years=log(sub_data$nb_years))

       

# Running the Stan model
mod <- sampling(stancode, data = stanlist, iter = 2000, chains = 4, cores = 4,init=0, save_warmup = FALSE)

# Save the model and the stanlist
list(mod = mod, stanlist = stanlist) %>% 
  saveRDS(file=here(paste0("outputs/ValidationNFI/StanModels/M62_",unique(sub_data$method_type),"_",unique(sub_data$method_input),".rds")))
    
# Model coefficients
conf95 <- broom.mixed::tidyMCMC(mod,
                                pars=params_to_estimate,
                                droppars = NULL, 
                                robust = FALSE, # give mean and standard deviation
                                ess = T, # effective sample size estimates
                                rhat = T, # Rhat estimates
                                conf.int = T, conf.level = 0.95 # include 95% credible intervals
                                ) %>% 
    dplyr::rename(conf95_low=conf.low,
                  conf95_high=conf.high,
                  mean=estimate,
                  std_deviation=std.error)
  
  
  broom.mixed::tidyMCMC(mod,
                        pars=params_to_estimate,
                        droppars = NULL, 
                        robust = TRUE, # give median and mean absolute deviation (= avg distance btw each data point and the mean)
                        ess = F, rhat = F, 
                        conf.int = T, conf.level = 0.80 # include 80% credible intervals
                        ) %>% 
    dplyr::rename(conf80_low=conf.low,
                  conf80_high=conf.high,
                  median=estimate,
                  mean_absolute_deviation=std.error) %>%
    inner_join(conf95,by=c("term")) %>% 
    mutate(method=method_i,
           method_type=unique(sub_data$method_type),
           method_input=unique(sub_data$method_input),
           .before=1)
    
  }) %>% bind_rows()

coefftab %>% saveRDS(file=here("outputs/ValidationNFI/ModelM62_coefftab.rds"))

Predictions of tree mortality probability

Plots

Let’s visualize the uncertainty around the estimation of the probability of tree mortality \(p\).

Code

############################################################################################
# Visualizing the relationships between GO predictions and mortality probability predictions
############################################################################################
df <- df %>% mutate(method_id=paste0(method_type,"_",method_input))

go_methods <- list(LFMM_all = "GO predictions from LFMM and all SNPs",
                LFMM_cand = "GO predictions from LFMM and candidate SNPs",
                LFMM_control = "GO predictions from LFMM and control SNPs",
                RDA_cand = "GO predictions from RDA and candidate SNPs",
                RDA_control = "GO predictions from RDA and control SNPs",
                GDM_cand = "GO predictions from GDM and candidate SNPs",
                GDM_control = "GO predictions from GDM and control SNPs",
                GF_cand = "GO predictions from GF and candidate SNPs",
                GF_control = "GO predictions from GF and control SNPs")
                
myplots <- lapply(unique(df$method_id), function(method_i){

# We load the stanlist
stanlist <- readRDS(file = here(paste0("outputs/ValidationNFI/StanModels/M62_",method_i,".rds")))$stanlist

# We load the model
mod <- readRDS(file = here(paste0("outputs/ValidationNFI/StanModels/M62_",method_i,".rds")))$mod

# We extract the posterior distributions of all parameters
post <- as.data.frame(mod)

# Vector of predictor values (based on the min and max of the GO predictions)
x_vec <- seq(min(stanlist$GO),max(stanlist$GO),0.05)

# Function to predict the mortality probability p with the initial tree height fixed to the mean
p_pred <- lapply(c("alpha_country[1]","alpha_country[2]"), function(country_i){
  
  f_p <- function(x) 1 -(exp(-exp(post[,country_i] + 
                                    post$`beta_DBH` * mean(stanlist$DBH) + 
                                    post$`beta_C` * mean(stanlist$C) + 
                                    post$`beta_C_DBH` * mean(stanlist$C) * mean(stanlist$DBH) + 
                                    post$`beta_GO` * x +
                                    log(5))))
  
sapply(x_vec, f_p) %>% 
    tibble::as_tibble() %>%
    rename_all(function(x) x_vec) %>%
    mutate(Iter = row_number()) %>%
    gather(x, p, -Iter) %>%
    group_by(x) %>%
    mutate(hpdi_l = HDInterval::hdi(p, credMass = 0.90)[1],
           hpdi_h = HDInterval::hdi(p, credMass = 0.90)[2],
           p_mean = mean(p)) %>%
    ungroup() %>%
    mutate(x = as.numeric(x),
           Country=ifelse(country_i=="alpha_country[1]","Spain","France"))
  
}) %>% bind_rows()


# Keep mean and 90% HDPI of p (one value for each iteration)
p_mean_df <- p_pred %>%
  dplyr::select(x,hpdi_l,hpdi_h,p_mean,Country) %>% 
  dplyr::distinct()


# Plots
#=======

# Plot options 
y_limits <- c(0,0.2)
x_axis <- "Mean-centered GO predictions"
y_axis <- "Probability of tree mortality"

p <- ggplot() +
  ylim(y_limits) +
  labs(y=y_axis, x=x_axis) +
  theme_bw() +
  theme(legend.position=c(0.8,0.85),
        legend.text = element_text(size=13),
        legend.title = element_text(size=13),
        legend.background = element_blank())

# First 100 posterior draws of the mortality probability p
p1 <- p +
  geom_point(data = p_pred %>% filter(Iter < 101),
             aes(x, p, color=Country), alpha = .2) 
  
# Mean (line) and 90% HPDI (shaded region) of the mortality probability p
p2 <- p +
  geom_ribbon(data = p_mean_df,
              aes(x = x, ymin = hpdi_l, ymax = hpdi_h, fill=Country),
              alpha = .2) +
  geom_line(data = p_mean_df,
              aes(x=x, y=p_mean,col=Country))

p1_p2 <- ggarrange(p1, p2 + labs(y=""), nrow = 1)

# Title
title <- ggdraw() + 
  draw_label(paste0(go_methods[[method_i]]),
             fontface = 'bold',
             x = 0, hjust = 0) +   
  theme(plot.margin = margin(0, 0, 0, 7))

# merge title and plots
p1_p2 <- plot_grid(title, p1_p2, ncol = 1, rel_heights = c(0.1, 1))

p2 <- p2 + labs(subtitle = go_methods[[method_i]])

ggsave(p1_p2,
       file=here(paste0("figs/ValidationNFI/MortalityProbabilityPredictions_UncertaintyIntervals/M62_",method_i,".pdf")),
       device="pdf",
       height=6,
       width=11)


return(p2)

})

p <- plot_grid(myplots[[1]] + labs(x="") + theme(legend.position = "none"),
               myplots[[2]] + labs(y="",x="") + theme(legend.position = "none"),
               myplots[[3]] + labs(y="",x=""),
               myplots[[4]] + labs(x="") + theme(legend.position = "none"),
               myplots[[5]] + labs(y="",x="") + theme(legend.position = "none"),
               myplots[[6]] + labs(y="",x="") + theme(legend.position = "none"),
               myplots[[7]] + theme(legend.position = "none"),
               myplots[[8]] + labs(y="") + theme(legend.position = "none"),
               myplots[[9]] + labs(y="") + theme(legend.position = "none"),
               ncol=3)

ggsave(p,
       file=here(paste0("figs/ValidationNFI/MortalityProbabilityPredictions_M62.pdf")),
       device="pdf",
       height=12,
       width=12)

Table

In the table below:

• mean and 90% credible intervals (highest posterior density intervals) of \(p\) for \(x=\{-1, 0, 1\}\).

• percent change in \(p\) associated with a one-unit increase in \(x\) (which corresponds to a one-standard deviation increase).

Code

p_pred <- lapply(unique(df$method_id), function(method_i){

# We load the stanlist
stanlist <- readRDS(file = here(paste0("outputs/ValidationNFI/StanModels/M62_",method_i,".rds")))$stanlist

# We load the model
mod <- readRDS(file = here(paste0("outputs/ValidationNFI/StanModels/M62_",method_i,".rds")))$mod

# We extract the posterior distributions of all parameters
post <- as.data.frame(mod)

# Vector of predictor values (based on the min and max of the GO predictions)
x_vec <- c(-1,0,1)

# Function to predict the mortality probability p with the initial tree height fixed to the mean
p_pred <- lapply(c("alpha_country[1]","alpha_country[2]"), function(country_i){
  
  f_p <- function(x) 1 -(exp(-exp(post[,country_i] + 
                                    post$`beta_DBH` * mean(stanlist$DBH) + 
                                    post$`beta_C` * mean(stanlist$C) + 
                                    post$`beta_C_DBH` * mean(stanlist$C) * mean(stanlist$DBH) + 
                                    post$`beta_GO` * x +
                                    log(5))))
  
sapply(x_vec, f_p) %>% 
    tibble::as_tibble() %>%
    rename_all(function(x) x_vec) %>%
    mutate(Iter = row_number()) %>%
    gather(x, p, -Iter) %>%
    group_by(x) %>%
    mutate(hpdi_l = HDInterval::hdi(p, credMass = 0.90)[1],
           hpdi_h = HDInterval::hdi(p, credMass = 0.90)[2],
           p_mean = mean(p)) %>%
    ungroup() %>%
    mutate(x = as.numeric(x),
           Country=ifelse(country_i=="alpha_country[1]","Spain","France"))
  
}) %>%  
  bind_rows() %>% 
  dplyr::select(x,p_mean,hpdi_l,hpdi_h,Country) %>% 
  distinct() %>% 
  mutate_if(is.numeric, round,  digits=3)

lapply(unique(p_pred$Country), function(country_i){
  
p_pred_sub <- p_pred %>% dplyr::filter(Country==country_i)

delta_p_pos <- (p_pred_sub$p_mean[p_pred_sub$x==1] * 100 / p_pred_sub$p_mean[p_pred_sub$x==0])-100
delta_p_neg <- (p_pred_sub$p_mean[p_pred_sub$x==-1] * 100 / p_pred_sub$p_mean[p_pred_sub$x==0])-100

tibble("Method" = go_methods[[method_i]],
       "Country"=country_i,
         "p(GO=-1)"=paste0(p_pred_sub$p_mean[p_pred_sub$x==-1], " [",
                       p_pred_sub$hpdi_l[p_pred_sub$x==-1],";",
                       p_pred_sub$hpdi_h[p_pred_sub$x==-1],"]"),
       "p(GO=0)"=paste0(p_pred_sub$p_mean[p_pred_sub$x==0], " [",
                       p_pred_sub$hpdi_l[p_pred_sub$x==0],";",
                       p_pred_sub$hpdi_h[p_pred_sub$x==0],"]"),
       "p(GO=1)"=paste0(p_pred_sub$p_mean[p_pred_sub$x==1], " [",
                       p_pred_sub$hpdi_l[p_pred_sub$x==1],";",
                       p_pred_sub$hpdi_h[p_pred_sub$x==1],"]"),
       "+Δp"=paste0(round(delta_p_pos,3),"%"),
       "-Δp"=paste0(round(delta_p_neg,3),"%"))
  
}) %>% bind_rows()

})  %>% bind_rows()

p_pred %>% saveRDS(file = here("tables/MortalityPredNFI.rds"))

Code

readRDS(file = here("tables/MortalityPredNFI.rds"))  %>% kable_mydf()

Method	Country	p(GO=-1)	p(GO=0)	p(GO=1)	+Δp	-Δp
GO predictions from GDM and control SNPs	Spain	0.129 [0.127;0.132]	0.116 [0.115;0.117]	0.104 [0.103;0.106]	-10.345%	11.207%
GO predictions from GDM and control SNPs	France	0.052 [0.049;0.054]	0.046 [0.044;0.048]	0.041 [0.039;0.043]	-10.87%	13.043%
GO predictions from GF and control SNPs	Spain	0.104 [0.102;0.105]	0.111 [0.11;0.112]	0.119 [0.117;0.12]	7.207%	-6.306%
GO predictions from GF and control SNPs	France	0.048 [0.046;0.05]	0.052 [0.05;0.054]	0.055 [0.053;0.058]	5.769%	-7.692%
GO predictions from RDA and control SNPs	Spain	0.117 [0.115;0.119]	0.113 [0.112;0.115]	0.109 [0.108;0.111]	-3.54%	3.54%
GO predictions from RDA and control SNPs	France	0.051 [0.049;0.054]	0.05 [0.047;0.052]	0.048 [0.046;0.05]	-4%	2%
GO predictions from LFMM and control SNPs	Spain	0.11 [0.108;0.111]	0.112 [0.111;0.113]	0.114 [0.113;0.116]	1.786%	-1.786%
GO predictions from LFMM and control SNPs	France	0.05 [0.047;0.052]	0.051 [0.048;0.053]	0.052 [0.049;0.054]	1.961%	-1.961%
GO predictions from GDM and candidate SNPs	Spain	0.087 [0.086;0.088]	0.108 [0.106;0.109]	0.133 [0.131;0.134]	23.148%	-19.444%
GO predictions from GDM and candidate SNPs	France	0.043 [0.041;0.044]	0.053 [0.05;0.055]	0.066 [0.063;0.069]	24.528%	-18.868%
GO predictions from GF and candidate SNPs	Spain	0.092 [0.091;0.094]	0.11 [0.109;0.111]	0.131 [0.13;0.133]	19.091%	-16.364%
GO predictions from GF and candidate SNPs	France	0.043 [0.041;0.045]	0.051 [0.049;0.053]	0.061 [0.059;0.064]	19.608%	-15.686%
GO predictions from RDA and candidate SNPs	Spain	0.099 [0.098;0.101]	0.11 [0.108;0.111]	0.121 [0.12;0.123]	10%	-10%
GO predictions from RDA and candidate SNPs	France	0.047 [0.045;0.049]	0.052 [0.05;0.055]	0.058 [0.056;0.061]	11.538%	-9.615%
GO predictions from LFMM and candidate SNPs	Spain	0.11 [0.108;0.111]	0.112 [0.111;0.113]	0.114 [0.113;0.116]	1.786%	-1.786%
GO predictions from LFMM and candidate SNPs	France	0.05 [0.047;0.052]	0.051 [0.049;0.053]	0.052 [0.05;0.055]	1.961%	-1.961%
GO predictions from LFMM and all SNPs	Spain	0.108 [0.106;0.11]	0.112 [0.11;0.113]	0.116 [0.114;0.117]	3.571%	-3.571%
GO predictions from LFMM and all SNPs	France	0.049 [0.047;0.051]	0.051 [0.049;0.054]	0.053 [0.051;0.056]	3.922%	-3.922%

Model coefficients

Below, we plot the 95% credible intervals of:

• the \(\beta_{GO}\) coefficients, which stand for the association between mortality rates in NFI plots and the genomic offset predictions (i.e. capturing the potential maladaptation of the populations at the location of the NFI plots).

• the \(\beta_C\) coefficients, which stands for the association between the tree basal area (all species confounded) in the NFI plots and the mortality rates.

  • the \(\beta_{DBH}\) coefficients, which stands for the association between the mean DBH (diameter at breast height) of all maritime pines in the NFI plots (including adults, dead trees and seedlings/saplings) and the mortality rates.

  • the \(\beta_{C*DBH}\) coefficients, which stands for the interaction between the mean DBH and the tree basal area.

Code

coefftab <- readRDS(file=here("outputs/ValidationNFI/ModelM62_coefftab.rds"))
  
coeff_match <- list(beta_GO=list("Regression coefficients $\\beta_{GO}$",c(0.77,0.85)),#"$\\beta_{GO}$ estimates",
                    beta_C=list("Regression coefficients $\\beta_C$",c(0.77,0.16)),
                    beta_DBH=list("Regression coefficients $\\beta_{DBH}$",c(0.77,0.16)),
                    beta_C_DBH=list("Regression coefficients $\\beta_{C*DBH}$",c(0.77,0.8)))

p <- lapply(c("beta_GO","beta_C","beta_DBH","beta_C_DBH"), function(coeff){
  
p <- coefftab %>% 
  filter(term==coeff) %>% 
  mutate(method_input=case_when(method_input=="all" ~ "All SNPs",
                                method_input=="cand" ~ "Candidate SNPs",
                                method_input=="control" ~ "Control SNPs")) %>% 
  mutate(method_input=factor(method_input, levels=c("All SNPs",
                                                    "Candidate SNPs",
                                                    "Control SNPs"))) %>% 
    ggplot(aes(x = method_type,
               y = median,
               ymin = conf95_low, 
               ymax = conf95_high,
               color=method_input,
               shape=method_input)) +
  {if(coeff=="beta_GO")  geom_rect(inherit.aes=FALSE, aes(xmin=-Inf, xmax=Inf, ymin=0, ymax=Inf), color="transparent", fill="green", alpha=0.01)} + 
  geom_hline(yintercept = 0,color="gray") +
  geom_pointinterval(position = position_dodge(width = .4),point_size=3.5,size=3) +
  ylab(TeX(coeff_match[[coeff]][[1]])) + xlab("") +
  ylim(c(-0.15,0.25)) +
  scale_color_manual(values=colors_coeff) +
  scale_shape_manual(values = c(rep(16,6),rep(17,3))) +
  theme_bw() +
  labs(color="SNP set",shape="SNP set") +
  theme(axis.text.x =  element_text(size=13),
        axis.text.y = element_text(size=13),
        axis.title.y = element_text(size=16),
        axis.title.x = element_text(size=1),
        legend.title=element_text(size=13), 
        legend.box.background = element_rect(colour = "gray", linewidth=0.6),
        strip.text.x = element_text(size = 16),
        strip.background = element_blank(),
        legend.position = coeff_match[[coeff]][[2]],
        legend.text=element_text(size=12),
        panel.grid.minor.x=element_blank(),
        panel.grid.major.x=element_blank()) +
  guides(color=guide_legend(ncol=1),
         shape = guide_legend(override.aes = list(size =2 )))

# save in pdf and png
ggsave(p,file=here(paste0("figs/ValidationNFI/ModelM62_",coeff,"_IntervalPlots.pdf")),
       device="pdf",
       height=5,
       width=5)

return(p)

})

Warning: Removed 1 rows containing missing values (`geom_segment()`).

Code

p

Warning: Removed 1 rows containing missing values (`geom_segment()`).

\(\beta_{GO}\) estimates against a null distribution

We evaluated the probability that the estimated relationship between the mortality rates and GO predictions could be obtained by chance. For that, we ran again the model on the NFI data but we replaced GO predictions by a randomly generated variable \(X\) such as \(X \sim \mathcal{N}(0,1)\). We repeated this step 100 times and obtained a distribution of \(\beta_X\) coefficients capturing the relationship between mortality rates and a randomly generated variable. We then compared this distribution of \(\beta_X\) coefficients to the estimated \(\beta_{GO}\) coefficients capturing the relationship between the mortality rates and the GO predictions.

Code

# Load the coefficient values of models with predicted GO
coefftab <- readRDS(file=here("outputs/ValidationNFI/ModelM62_coefftab.rds")) %>% 
  mutate(run_ID=method,
         GO="predicted") %>% 
  dplyr::filter(term == "beta_GO") %>%
  dplyr::select(run_ID, GO, mean, contains("conf"))

# Load the coefficient values of models with random GO and merge with the predicted GO coefficient table
coefftab <- readRDS(file=here(paste0("outputs/ValidationNFI/Simulations/m6_2_truedata_randomGO_100simulations.rds"))) %>% 
  bind_rows(.id="sim_ID") %>% 
  dplyr::filter(term == "beta_GO") %>%
  mutate(run_ID = paste0("sim ",sim_ID),#run_ID = paste0("Simulation ",sim_ID),
         GO="random") %>% 
  dplyr::select(run_ID, GO, mean, contains("conf")) %>% 
  bind_rows(coefftab)
  
# Plot the mean and 95% credible intervals with interval plots
p <- coefftab %>% 
  ggplot(aes(mean, reorder(run_ID, mean))) +
  geom_vline(xintercept = 0, color="gray30") +
  geom_errorbar(aes(xmin = conf95_low, xmax = conf95_high, color=GO)) +
  scale_color_manual(values=c("magenta","gray"),
                     name = "Genomic offset") +
  geom_point() +
  theme_bw() +
  labs(x="GO effect size", y="") +
  theme(axis.text.y = element_text(size=6))

# save the plot in pdf
ggsave(p,
       file=here(paste0("figs/ValidationNFI/ModelM62_IntervalPlots_RandomGO_vs_PredictedGO.pdf")),
       device="pdf",
       height=10,
       width=8)

p

GO predictions in Galicia

GO predictions in Galicia (especially those obtained with GF) show a strong boundary between two zones (see Figure 4b of the manuscript). How can we explain this pattern?

Code

nfi_data %>% 
  mutate(year_second_survey = ifelse(second_survey > 1999, "After 2000",second_survey)) %>% 
           ggplot() + 
    geom_sf(data = world, fill="gray98") + 
    theme_bw() +
    scale_x_continuous(limits = c(-10, 12)) +
    scale_y_continuous(limits = c(35, 51)) + 
    geom_point( aes(x=longitude,
                    y=latitude,
                    color=factor(year_second_survey)),
               alpha=0.4, size=0.5) + 
    xlab("") + ylab("") +
    ggtitle("") +
    scale_color_manual(name = "Year of the second survey", 
                       values =  c("yellow2","darkorchid","#60F2E7","#C1E5AC")) +
    theme(legend.box.background = element_rect(colour = "gray80", linewidth=0.6)) +
   guides(colour = guide_legend(override.aes = list(size=2)))

In Galicia, GO predictions from the GF approach (see Figure 4b of the manuscript) are higher for plots surveyed for the second time in 1998 than plots surveyed for the second time in 1997.

We explore the climatic differences between the two sets of plots, ie plots surveyed in 1997 vs 1998. First, we extract the annual climatic values at the location of the Galician plots sampled in 1997 and 1998 and plot their values across the study period (ie between the first and second surveys).

Code

# Extracting annual climatic values at the location of the NFI plots in Galicia surveyed for the second time inm 1997 or 1998
subclim <- clim %>% 
  left_join(nfi_data, by=c("plotcode","longitude","latitude")) %>% 
    dplyr::filter(second_survey %in% 1997:1998) %>% 
    dplyr::filter(year %in% min(first_survey):1998) %>% 
    dplyr::select(second_survey,year,all_of(clim_var))


# violin_plots <- lapply(clim_var, function(x) {
#   
# var_name <- extract_climatedt_metadata(var_clim = x) %>% 
#   mutate(var_legend= paste0(description, " (", label,"; ",unit_symbol,")")) %>% 
#   pull(var_legend)
# 
# subclim %>% 
#   mutate(year=as.factor(year)) %>% 
#   ggplot(aes_string(y=x,x="year")) + 
#   geom_jitter(shape=16,aes(colour=factor(second_survey)),position=position_jitter(0.2),size=2,alpha=0.3) +
#   xlab("") + 
#   scale_color_manual(name="Year of the second survey",
#                      values = c("yellow2","darkorchid")) +
#   ggtitle(var_name) +
#   ylab("mean climate under survey period - mean past climate") +
#   theme_bw()
# 
# })

violin_plots <- lapply(clim_var, function(x) {
  
var_name <- extract_climatedt_metadata(var_clim = x) %>% 
  mutate(var_legend= paste0(description, " (", label,"; ",unit_symbol,")")) %>% 
  pull(var_legend)

subclim %>% 
  mutate(year=as.factor(year),
         second_survey=as.factor(second_survey)) %>% 
  ggplot(aes_string(y=x,x="year",color="second_survey")) + 
  geom_point(shape=16,position=position_jitter(0.2),size=2,alpha=0.1) +
  geom_violin(alpha=0.3,fill="white", linewidth=0.7) + 
  xlab("") + ylab("") +
  ggtitle(var_name) +
  scale_color_manual(name="Year of the second survey",
                     values = c("yellow2","darkorchid")) +
  theme_bw() +
  guides(colour = guide_legend(override.aes = list(size=2)))

})

Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
ℹ Please use tidy evaluation idioms with `aes()`.
ℹ See also `vignette("ggplot2-in-packages")` for more information.

Code

violin_plots 

From these plots, it seems that the plots sampled in 1997 have lower isothermality (bio3) and lower temperature seasonality (bio4).

We also compare the mean climatic values of the reference and survey periods between the plots surveyed in 1997 and 1998. For that, we calculate for each plot (and each climatic variable) the difference between its mean climate across the reference period (ie 1901-1950) and its mean climate across the survey period (between the first and second survey period).

Code

# We load the datasets with mean climatic values across the reference and survey periods
# Then we calculate the difference between the mean climates of the two periods for each plot
#  We then keep only plots surveyed (second survey) in 1997 or 1998 (Galician plots)
subnficlim <- readRDS(here("data/ClimaticData/NFIplots/NFIclimate.rds")) %>% 
  bind_rows() %>% 
  group_by(plotcode,longitude,latitude,elevation) %>% 
  summarise(across(tmn01:Eref, diff)) %>% 
  ungroup() %>% 
  left_join(nfi_data, by=c("plotcode","longitude","latitude")) %>%
  dplyr::filter(second_survey %in% 1997:1998)

bg_color <- "gray44" # background color

violin_plots <- lapply(clim_var, function(x) {
  
var_name <- extract_climatedt_metadata(var_clim = x) %>% 
  mutate(var_legend= paste0(description, " (", label,"; ",unit_symbol,")")) %>% 
  pull(var_legend)

subnficlim %>% 
  mutate(second_survey=as.factor(second_survey)) %>% 
  ggplot(aes_string(y=x,x="second_survey")) + 
  geom_violin(alpha=0.2, fill="gray76", color="darkgreen") + 
  geom_jitter(shape=16,position=position_jitter(0.2),size=2,alpha=0.4) +
  xlab("") + 
  ggtitle(var_name) +
  ylab("mean climate under survey period - mean past climate") +
  theme_bw()

})

violin_plots 

These plots show that plots surveyed in 1998 show higher increase in the mean annual temperature (bio1), lower increase in annual precipitation (bio12) and stronger decrease in temperature seasonality (bio4) than plots surveyed in 1997. As temperature seasonality is the most important variable to explain the allelic turnover in GF, it may explain the strong differences in GO predictions among plots surveyed in 1997 and 1998.

Session information

Code

devtools::session_info()

─ Session info ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
 setting  value
 version  R version 4.3.2 (2023-10-31)
 os       Ubuntu 22.04.3 LTS
 system   x86_64, linux-gnu
 ui       X11
 language (EN)
 collate  fr_FR.UTF-8
 ctype    fr_FR.UTF-8
 tz       Europe/Paris
 date     2024-02-01
 pandoc   3.1.1 @ /usr/lib/rstudio/resources/app/bin/quarto/bin/tools/ (via rmarkdown)

─ Packages ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
 package           * version  date (UTC) lib source
 abind               1.4-5    2016-07-21 [1] CRAN (R 4.3.2)
 arrayhelpers        1.1-0    2020-02-04 [1] CRAN (R 4.3.2)
 backports           1.4.1    2021-12-13 [1] CRAN (R 4.3.2)
 bit                 4.0.5    2022-11-15 [1] CRAN (R 4.3.2)
 bit64               4.0.5    2020-08-30 [1] CRAN (R 4.3.2)
 broom               1.0.5    2023-06-09 [1] CRAN (R 4.3.2)
 cachem              1.0.8    2023-05-01 [1] CRAN (R 4.3.2)
 car                 3.1-2    2023-03-30 [1] CRAN (R 4.3.2)
 carData             3.0-5    2022-01-06 [1] CRAN (R 4.3.2)
 cellranger          1.1.0    2016-07-27 [1] CRAN (R 4.3.2)
 checkmate           2.3.1    2023-12-04 [1] CRAN (R 4.3.2)
 class               7.3-22   2023-05-03 [4] CRAN (R 4.3.1)
 classInt            0.4-10   2023-09-05 [1] CRAN (R 4.3.2)
 cli                 3.6.2    2023-12-11 [1] CRAN (R 4.3.2)
 coda                0.19-4   2020-09-30 [1] CRAN (R 4.3.2)
 codetools           0.2-19   2023-02-01 [4] CRAN (R 4.2.2)
 colorspace          2.1-0    2023-01-23 [1] CRAN (R 4.3.2)
 cowplot           * 1.1.2    2023-12-15 [1] CRAN (R 4.3.2)
 crayon              1.5.2    2022-09-29 [1] CRAN (R 4.3.2)
 DBI                 1.2.0    2023-12-21 [1] CRAN (R 4.3.2)
 devtools            2.4.5    2022-10-11 [1] CRAN (R 4.3.2)
 digest              0.6.34   2024-01-11 [1] CRAN (R 4.3.2)
 distributional      0.3.2    2023-03-22 [1] CRAN (R 4.3.2)
 dplyr             * 1.1.4    2023-11-17 [1] CRAN (R 4.3.2)
 e1071               1.7-14   2023-12-06 [1] CRAN (R 4.3.2)
 ellipsis            0.3.2    2021-04-29 [1] CRAN (R 4.3.2)
 evaluate            0.23     2023-11-01 [1] CRAN (R 4.3.2)
 fansi               1.0.6    2023-12-08 [1] CRAN (R 4.3.2)
 farver              2.1.1    2022-07-06 [1] CRAN (R 4.3.2)
 fastmap             1.1.1    2023-02-24 [1] CRAN (R 4.3.2)
 forcats           * 1.0.0    2023-01-29 [1] CRAN (R 4.3.2)
 fs                  1.6.3    2023-07-20 [1] CRAN (R 4.3.2)
 generics            0.1.3    2022-07-05 [1] CRAN (R 4.3.2)
 geodist           * 0.0.8    2022-10-19 [1] CRAN (R 4.3.2)
 ggdist              3.3.1    2023-11-27 [1] CRAN (R 4.3.2)
 ggplot2           * 3.4.4    2023-10-12 [1] CRAN (R 4.3.2)
 ggpubr            * 0.6.0    2023-02-10 [1] CRAN (R 4.3.2)
 ggrepel           * 0.9.4    2023-10-13 [1] CRAN (R 4.3.2)
 ggsignif            0.6.4    2022-10-13 [1] CRAN (R 4.3.2)
 glue                1.6.2    2022-02-24 [1] CRAN (R 4.3.2)
 goftest             1.2-3    2021-10-07 [1] CRAN (R 4.3.2)
 gridExtra           2.3      2017-09-09 [1] CRAN (R 4.3.2)
 gtable              0.3.4    2023-08-21 [1] CRAN (R 4.3.2)
 HDInterval          0.2.4    2022-11-17 [1] CRAN (R 4.3.2)
 here              * 1.0.1    2020-12-13 [1] CRAN (R 4.3.2)
 highr               0.10     2022-12-22 [1] CRAN (R 4.3.2)
 hms                 1.1.3    2023-03-21 [1] CRAN (R 4.3.2)
 htmltools           0.5.7    2023-11-03 [1] CRAN (R 4.3.2)
 htmlwidgets         1.6.4    2023-12-06 [1] CRAN (R 4.3.2)
 httpuv              1.6.13   2023-12-06 [1] CRAN (R 4.3.2)
 httr                1.4.7    2023-08-15 [1] CRAN (R 4.3.2)
 inline              0.3.19   2021-05-31 [1] CRAN (R 4.3.2)
 janitor           * 2.2.0    2023-02-02 [1] CRAN (R 4.3.2)
 jsonlite            1.8.8    2023-12-04 [1] CRAN (R 4.3.2)
 kableExtra        * 1.3.4    2021-02-20 [1] CRAN (R 4.3.2)
 KernSmooth          2.23-22  2023-07-10 [4] CRAN (R 4.3.1)
 knitr             * 1.45     2023-10-30 [1] CRAN (R 4.3.2)
 labeling            0.4.3    2023-08-29 [1] CRAN (R 4.3.2)
 later               1.3.2    2023-12-06 [1] CRAN (R 4.3.2)
 latex2exp         * 0.9.6    2022-11-28 [1] CRAN (R 4.3.2)
 lattice             0.22-5   2023-10-24 [4] CRAN (R 4.3.1)
 lifecycle           1.0.4    2023-11-07 [1] CRAN (R 4.3.2)
 lmom                3.0      2023-08-29 [1] CRAN (R 4.3.2)
 lmomco              2.4.11   2023-08-30 [1] CRAN (R 4.3.2)
 Lmoments            1.3-1    2019-03-15 [1] CRAN (R 4.3.2)
 loo                 2.6.0    2023-03-31 [1] CRAN (R 4.3.2)
 lubridate         * 1.9.3    2023-09-27 [1] CRAN (R 4.3.2)
 magrittr          * 2.0.3    2022-03-30 [1] CRAN (R 4.3.2)
 MASS                7.3-60   2023-05-04 [4] CRAN (R 4.3.1)
 matrixStats         1.2.0    2023-12-11 [1] CRAN (R 4.3.2)
 memoise             2.0.1    2021-11-26 [1] CRAN (R 4.3.2)
 mime                0.12     2021-09-28 [1] CRAN (R 4.3.2)
 miniUI              0.1.1.1  2018-05-18 [1] CRAN (R 4.3.2)
 munsell             0.5.0    2018-06-12 [1] CRAN (R 4.3.2)
 paletteer         * 1.5.0    2022-10-19 [1] CRAN (R 4.3.2)
 pillar              1.9.0    2023-03-22 [1] CRAN (R 4.3.2)
 pkgbuild            1.4.3    2023-12-10 [1] CRAN (R 4.3.2)
 pkgconfig           2.0.3    2019-09-22 [1] CRAN (R 4.3.2)
 pkgload             1.3.3    2023-09-22 [1] CRAN (R 4.3.2)
 plyr                1.8.9    2023-10-02 [1] CRAN (R 4.3.2)
 posterior           1.5.0    2023-10-31 [1] CRAN (R 4.3.2)
 profvis             0.3.8    2023-05-02 [1] CRAN (R 4.3.2)
 promises            1.2.1    2023-08-10 [1] CRAN (R 4.3.2)
 proxy               0.4-27   2022-06-09 [1] CRAN (R 4.3.2)
 purrr             * 1.0.2    2023-08-10 [1] CRAN (R 4.3.2)
 QuickJSR            1.0.9    2023-12-18 [1] CRAN (R 4.3.2)
 R6                  2.5.1    2021-08-19 [1] CRAN (R 4.3.2)
 ragg                1.2.7    2023-12-11 [1] CRAN (R 4.3.2)
 raster              3.6-26   2023-10-14 [1] CRAN (R 4.3.2)
 Rcpp                1.0.11   2023-07-06 [1] CRAN (R 4.3.2)
 RcppParallel        5.1.7    2023-02-27 [1] CRAN (R 4.3.2)
 readr             * 2.1.4    2023-02-10 [1] CRAN (R 4.3.2)
 readxl            * 1.4.3    2023-07-06 [1] CRAN (R 4.3.2)
 rematch2            2.1.2    2020-05-01 [1] CRAN (R 4.3.2)
 remotes             2.4.2.1  2023-07-18 [1] CRAN (R 4.3.2)
 reshape             0.8.9    2022-04-12 [1] CRAN (R 4.3.2)
 reshape2          * 1.4.4    2020-04-09 [1] CRAN (R 4.3.2)
 rlang               1.1.3    2024-01-10 [1] CRAN (R 4.3.2)
 rmarkdown           2.25     2023-09-18 [1] CRAN (R 4.3.2)
 rnaturalearth     * 1.0.1    2023-12-15 [1] CRAN (R 4.3.2)
 rnaturalearthdata   0.1.0    2017-02-21 [1] CRAN (R 4.3.2)
 rprojroot           2.0.4    2023-11-05 [1] CRAN (R 4.3.2)
 rstan             * 2.32.3   2023-10-15 [1] CRAN (R 4.3.2)
 rstatix             0.7.2    2023-02-01 [1] CRAN (R 4.3.2)
 rstudioapi          0.15.0   2023-07-07 [1] CRAN (R 4.3.2)
 rvest               1.0.3    2022-08-19 [1] CRAN (R 4.3.2)
 scales              1.3.0    2023-11-28 [1] CRAN (R 4.3.2)
 sessioninfo         1.2.2    2021-12-06 [1] CRAN (R 4.3.2)
 sf                  1.0-15   2023-12-18 [1] CRAN (R 4.3.2)
 shiny               1.8.0    2023-11-17 [1] CRAN (R 4.3.2)
 snakecase           0.11.1   2023-08-27 [1] CRAN (R 4.3.2)
 sp                  2.1-2    2023-11-26 [1] CRAN (R 4.3.2)
 SPEI              * 1.8.1    2023-03-02 [1] CRAN (R 4.3.2)
 StanHeaders       * 2.26.28  2023-09-07 [1] CRAN (R 4.3.2)
 stringi             1.8.3    2023-12-11 [1] CRAN (R 4.3.2)
 stringr           * 1.5.1    2023-11-14 [1] CRAN (R 4.3.2)
 svglite             2.1.3    2023-12-08 [1] CRAN (R 4.3.2)
 svUnit              1.0.6    2021-04-19 [1] CRAN (R 4.3.2)
 systemfonts         1.0.5    2023-10-09 [1] CRAN (R 4.3.2)
 tensorA             0.36.2.1 2023-12-13 [1] CRAN (R 4.3.2)
 terra               1.7-65   2023-12-15 [1] CRAN (R 4.3.2)
 textshaping         0.3.7    2023-10-09 [1] CRAN (R 4.3.2)
 tibble            * 3.2.1    2023-03-20 [1] CRAN (R 4.3.2)
 tidybayes         * 3.0.6    2023-08-12 [1] CRAN (R 4.3.2)
 tidyr             * 1.3.0    2023-01-24 [1] CRAN (R 4.3.2)
 tidyselect          1.2.0    2022-10-10 [1] CRAN (R 4.3.2)
 tidyverse         * 2.0.0    2023-02-22 [1] CRAN (R 4.3.2)
 timechange          0.2.0    2023-01-11 [1] CRAN (R 4.3.2)
 TLMoments           0.7.5.3  2022-03-27 [1] CRAN (R 4.3.2)
 tzdb                0.4.0    2023-05-12 [1] CRAN (R 4.3.2)
 units               0.8-5    2023-11-28 [1] CRAN (R 4.3.2)
 urlchecker          1.0.1    2021-11-30 [1] CRAN (R 4.3.2)
 usethis             2.2.2    2023-07-06 [1] CRAN (R 4.3.2)
 utf8                1.2.4    2023-10-22 [1] CRAN (R 4.3.2)
 vctrs               0.6.5    2023-12-01 [1] CRAN (R 4.3.2)
 viridisLite         0.4.2    2023-05-02 [1] CRAN (R 4.3.2)
 vroom               1.6.5    2023-12-05 [1] CRAN (R 4.3.2)
 webshot             0.5.5    2023-06-26 [1] CRAN (R 4.3.2)
 withr               3.0.0    2024-01-16 [1] CRAN (R 4.3.2)
 xfun                0.41     2023-11-01 [1] CRAN (R 4.3.2)
 xml2                1.3.6    2023-12-04 [1] CRAN (R 4.3.2)
 xtable            * 1.8-4    2019-04-21 [1] CRAN (R 4.3.2)
 yaml                2.3.8    2023-12-11 [1] CRAN (R 4.3.2)
 zoo                 1.8-12   2023-04-13 [1] CRAN (R 4.3.2)

 [1] /home/juliette/R/x86_64-pc-linux-gnu-library/4.3
 [2] /usr/local/lib/R/site-library
 [3] /usr/lib/R/site-library
 [4] /usr/lib/R/library

──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────

References

Changenet, Alexandre, Paloma Ruiz-Benito, Sophia Ratcliffe, Thibaut Fréjaville, Juliette Archambeau, Annabel J Porte, Miguel A Zavala, Jonas Dahlgren, Aleksi Lehtonen, and Marta Benito Garzón. 2021. “Occurrence but Not Intensity of Mortality Rises Towards the Climatic Trailing Edge of Tree Species Ranges in European Forests.” Global Ecology and Biogeography.